AEI-2004-050 



The Hough transform search for continuous gravitational waves 

Badri Krishnan,^'* Alicia M. Sintes,^' ^' ^ Maria Alessandra Papa,^'-!- 
Bernard F. Schutz,^' § Sergio Frasca,''' ^' ^ and Cristiano Palomba^^ ** 

^ Max-Planck-Institut fur Gravitationsphysik, Albert- Einstein- Instttut, Am Miihlenberg 1, D-14-476 Golm, Germany 
Departament de Fisica, Universitat de les Illes Balears, Cra. Valldemossa Km. 7.5, E-07122 Palma de Mallorca, Spain 

'^Cardiff University, Cardiff, CF2 3YB, United Kingdom 
^Dip. di Fisica, Universita di Roma "La Sapienza" P.le A. Mora 2, 1-00185 Rome, Italy 
^INFN Sezione di Roma, P.le A. Moro 2, 1-00185 Rome, Italy 

This paper describes an incoherent method to search for continuous gravitational waves based 
on the Hough transform, a well known technique used for detecting patterns in digital images. We 
apply the Hough transform to detect patterns in the time-frequency plane of the data produced 
by an earth-based gravitational wave detector. Two different flavors of searches will be considered, 
depending on the type of input to the Hough transform: either Fourier transforms of the detector 
data or the output of a coherent matched-filtering type search. We present the technical details for 
implementing the Hough transform algorithm for both kinds of searches, their statistical properties, 
and their sensitivities. 
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I. INTRODUCTION 

Rapidly rotating neutron stars are expected to be the 
primary sources of continuous gravitational waves, and 
the current generation of earth-based gravitational wave 
detectors might be able to detect them. Recent analysis 
of data from the first science runs of the LIGO [1-3] and 
GEO [1, 4, 5] interferometric detectors has already led 
to upper limits on the gravitational waves emitted by 
the pulsar J1939-I-2134 and its equatorial ellipticity [6]. 
The analysis of future science runs is expected to lead to 
upper limits below other astrophysical constraints, and 
eventually to detections. 

The analyses presented in [6] were based on the 
coherent integration of the detectors' output for the 
entire observation time (approximately 17 days) and 
used a Bayesian time-domain method and a frequentist 
frequency-domain [7] approach. The searches were not 
computationally expensive, targeting a single known pul- 
sar and processing only a narrow frequency band of about 
0.5 Hz around the pulsar emission frequency for a fixed 
sky location and spin-down rate known from radio obser- 
vations. 

Future continuous wave searches will involve search- 
ing longer data stretches (of order weeks to months) for 
unknown sources over a large frequency band, vast por- 
tions of the sky and spin-down parameter values. It 
is well known that the computational cost of coherent 
techniques for searches of this type is absolutely pro- 
hibitive [8]. Thus hierarchical methods have been pro- 
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posed. 

In hierarchical strategies incoherent techniques (less 
sensitive and less computationally expensive) are used 
to scan the data and the parameter space for interest- 
ing candidates which are then followed up with coherent 
searches. Different strategies can be envisaged that com- 
bine the data incoherently. All methods use, in some way, 
the power from the Fourier transforms of short stretches 
of data: in the frequency bins where the signal is present 
there should systematically be an excess of power. In 
order to compensate for the frequency modulation im- 
posed on the signal by the Earth's motion and the pul- 
sar's spin-down during the observation period, one must 
use not the power from the same frequency bins in each 
successive Fourier transform, but rather from the bins 
where one expects the signal peak to be. 

In the so called stack-slide method, one "slides" the 
frequency bins of each Fourier transform to line-up the 
signal peaks and then simply sums the power [9]. The 
Hough transform method can be seen as a variation on 
this where, after the sliding, one sums not the power but 
just zeros and ones, depending on whether the power 
in the frequency bin exceeds a threshold or meets some 
other criterion. Whereas in low signal-to-noise conditions 
in Gaussian noise, the standard power summing method 
is possibly optimal, the Hough transform method might 
be more robust in the presence of large spectral distur- 
bances. To see this, consider the case when a large spec- 
tral disturbance is present only in a single Fourier trans- 
form. This could have a very large effect on the power 
sum statistic, but no matter how large, this spectral line 
could only add -f 1 to the Hough statistic. 

The Hough transform is a robust parameter estima- 
tor of multi-dimensional patterns in images and it finds 
many applications in astronomical data analysis [10-12]. 
In the context of image processing, it provides robustness 
against missing data points or discontinuous features [13]. 
It was initially developed by Paul Hough to analyze bub- 
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FIG. 1: The detector frame and the wave frame. 



ble chamber pictures at CERN, and later patented by 
IBM [14, 15]. It is currently being used to analyze data 
from the LIGO and GEO detectors. The codes employed 
for these analyses are freely available as part of the LIGO 
Algorithms Library [16]. The VIRGO project [17, 18] 
is also setting up a similar hierarchical search pipeline. 
Studies of hierarchical strategies can be found in [19-24]. 

This paper is organized as follows: section II briefly 
describes the expected waveforms from an isolated spin- 
ning neutron star and summarizes the general strategy 
of a hierarc;liic;al search. Section III presents the general 
idea of the Hough transform and section IV describes 
its implementation for non-demodulated input data, and 
section V studies its statistical properties. Section VI de- 
scribes the Hough search using demodulated input data 
and finally section VII summarizes our main results. 



II. PRELIMINARIES 
A. The signal from a pulsar 

In this subsection we fix our notation and briefly review 
the expected gravitational wave signal from a spinning 
neutron star. Further details about the pulsar signal can 
be found in [7]; a concise review of the possible physical 
mechanisms that may be causing pulsars to emit grav- 
itational waves can be found in [6]. For our purposes, 
we only need the form of the gravitational wave signal as 
seen by an Earth based detector. 

Let ni and n2 denote the unit vectors pointing along 
the arms of the detector and denote by ( the angle be- 
tween the arms. Let z be the unit vector parallel to 
ni X n2. Apart from the detector frame (ni,n2,z), we 
also have the wave frame (x^ , , ) in which the unit 
vector z^ is along the direction of propagation of the 
wave and (xu,, y^,, z^,) form a right-handed orthonormal 
system. Finally, n = — z^ is the unit vector pointing 
in the direction of the neutron star; see figure 1. The 
spacetime metric (/^^ can be written as a perturbation 
of the flat metric 7]^,^: g^^, = rj^i, + h^^. The received 



gravitational wave h^^ has the form 

= h+{t) (e+)^, + /ix(i) (ex)^, (2.1) 

where e+ = x^, (g) x„ - y^, ® and ex = x„ (g) y^ -|- 
Yw ® x^, and t denotes clock time at the location of 
the (moving, accelerating) detector, which we refer to as 
detector time. The waveforms for the two polarizations 
are 

h+{t) = A+ cos ^{t) , hx{t) = A^sm^{t) (2.2) 

where $(t) is the phase of the gravitational wave and 
A_|_^x are the amplitudes; A+,x are constant in time and 
depend on the other pulsar parameters such as its rota- 
tional frequency, moments of inertia, the orientation of 
its rotation axis, its distance from Earth etc. The phase 
^{t) takes its simplest form when the time coordinate 
used is ^NS) the proper time in the rest frame of the neu- 
tron star: 

$Ns(tKs) =ct>o + 2^E T;^*-' (2.3) 

n=0 ^ 

where (j)o, f^^^^^ and f^^^'' (n > 1) arc respectively the 
phase, instantaneous frequency and the spin-down pa- 
rameters in the rest frame of the star at the fiducial start 
time t^s = 0, and s is the number of spin-down parame- 
ters included in our search. 

We refer the reader to [7] for the expression of $(t) in 
the detector frame as a hmction of detector time. For our 
purposes, we only need to know that the instantaneous 
frequency f{t) of the wave as observed by the detector 
is given, to a very good approximation, by the familiar 
non-relativistic Doppler formula: 

where v{t) is the detector velocity in the Solar System 
Barycenter (SSB) frame and / is given by 

/w = /(« + t^('-'o + ^)" (") 

n=l ^ ' 

where to is the fiducial detector time at the start of the 
observation, the /(„) are the spin-down parameters as 
measured in the SSB frame (these need not be equal to 
the fl^\ see [7]), and Ar{t) := r{t) - r(to) with r{t) 
being the position of the detector in the SSB frame at 
time t. We have also assumed the neutron star to be 
moving with uniform speed relative to the Sun and is 
so far away that there are no observable proper-motion 
effects. (These could be taken into account if necessary, 
at the cost of introducing further parameters.) 

The detector output is a linear combination of h+ and 
hx- 

h{t) = F+(n, V) h+{t) + Fx (n, ^) (t) (2.6) 
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where i^+,x are known as the the antenna pattern func- 
tions of the detector and depend on the direction n to 
the star and also on the polarization angle ij) which deter- 
mines the orientation of the (x^, y^) axes in their plane. 
In addition, the antenna pattern functions also depend 
on the detector parameters such as its latitude, the angle 
C between its arms, and the azimuth of the bisector of 
the arms. Due to the motion of the Earth, F+. x (n, ip) de- 
pend implicitly on time and for notational convenience, 
we shall usually denote the antenna pattern functions as 
F+^x(i)- Thus, the received signal is both amplitude- 
and frequency-modulated. 

The search method described in this paper depends on 
finding a signal whose frequency evolution fits the pattern 
produced by the Doppler shift and the spin-down. The 
parameters which determine this pattern are the ones 
which appear in equation (2.4), namely, (/(o), {/(„)}, n); 

these parameters will be collectively denoted by ^. 

The amplitudes A+^x are determined by the other pul- 
sar parameters such as the orientation of its axis, its el- 
lipticity, its distance from Earth etc. The search method 
presented in this paper depends only on the phase model 
of equation (2.3). The exact form of the amplitudes is 
model dependent. As an illustrative example, consider 
the wave emitted by a deformed spinning neutron star as 
in [6]. If fr is the rotational frequency of the star, the fre- 
quency of the gravitational wave is 2/^. The additional 
parameters determining this component of the pulsar sig- 
nal are l and h{) where l is the angle between the pulsar's 
axis of rotation and the vector = — n, and charac- 
terizes the amplitude of the emitted gravitational wave. 
The amplitudes A+^x are: 



A+ = -/io(l + cos l) , 
j4x = /iQCOsi-. 



(2.7) 
(2.8) 



If we assume the emission mechanism is due to deviations 
of the pulsar's shape from perfect axial symmetry, then 
the amplitude /iq will be 



ho = 



c4 



d 



(2.9) 



where d is the distance of the star from Earth, I^z is 
the z-z component of the star's moment of inertia with 
the ^-axis being its spin axis, and e := {Ixx — Iyy)/Izz 
is the equatorial ellipticity of the star. Among all the 
quantities appearing in this equation, the value of e is by 
far the most uncertain. Typical values are expected to be 
~ 10~* for standard neutron stars and values of ~ 10~^ 
are expected to be the maximum values [8]. There is 
also a very small uncertainty in the value of fr because 
the pulsar could have a (presently unobservable) radial 
velocity. This would produce a Doppler shift between the 
true value of fr in the neutron star (NS) frame, and its 
measured value on Earth. Assuming typical values of the 
radial velocity to be the same as the typically measured 
transverse velocities (~ 500 km/s, see e. g. [25]), we get 
an uncertainty in fr of ~ 0.1%. 



B. A multi-stage hierarchical search 

Consider performing a blind search for pulsars using a 

bank of templates and relying only on coherent matched 
filter techniques. Since a larger observation time im- 
plies better resolution in the space of frequency, spin- 
downs and sky-positions, the number of templates in- 
creases rapidly as a function of the total observation 
time. A typical example is an all-sky search for young, 
fast pulsars, i.e. for hypothetical signals with frequency 
/ < /max = lOOOHz and spin-down ages greater than 
T > T„i„ — 40yr. Let s be the number of spin-down pa- 
rameters that we search over and let T^bs be the total 
observation time. The number of templates required for 
this search has been calculated in equation (6.3) of [8]: 



Np « max,g{o,i...}[^^s^s(7;bs)] 



where 



iV,, 



(f^ 
llkHz 



s+2 



40yr 



s(s-|-l)/2 



(2.10) 



(2.11) 



gives the spin-down scaling and Fg is a function that 
depends on the observation time; for large observation 
times, Fs oc T^^^^. We have taken the maximum allowed 
fractional mismatch in observed signal power between the 
signal and the template to be 0.3. For example, if ,s = 2, 
assuming the observation time is significantly longer than 
a day, equation (6.7) of [8] approximates to : 



F2{T,^s) « 2.2 X 10^ X 



\ldayj 



(2.12) 



Thus, even for a 10 day search over two spin-down param- 
eters, Np « 2 X 10^^. The computational requirements 
for a search over these many templates is also estimated 
in [8]. It turns out that for the 10 day long search, if we 
wished to analyze the data in roughly real time, we would 
require a computational power of ~ 10® GFlops; for ref- 
erence, the fastest supercomputers ca. 2004 can do 'only' 
~ 10'' GFlops. Even if we insisted on searching over only 
a single spin-down parameter, for an observation period 
of only 10 days, the computational requirement turns out 
to be ~ 10^ GFlops. We therefore conclude that a search 
over any significant portion of parameter space for un- 
known pulsars is not possible in the foreseeable future if 
we restrict ourselves to fully coherent methods. 

One possible way to perform such a blind search would 
be to exploit the fact that Np increases faster than lin- 
early with T„bs. Thus if we break up the data set into 
smaller segments, it might be feasible to analyze each 
data segment coherently. An incoherent method is then 
used as a computationally inexpensive and sub-optimal 
way of combining the outputs of the different coherent 
segments. This would be one step in a multi-stage hier- 
archical scheme; see figure 2. 

In this scheme, we start with a data stream covering 
a total observation time T^b^. Divide the available data 
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Acquire data 
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Break up data into 
smaller segments 



Analyze each segment 
coherently 



Acquire more data 



Combine segments 
incoherently 



Select candidates in 
parameter space 



Analyse candidates 
fully coherently 



Announce detection 
or set upper-limits 



FIG. 2: A hierarchical scheme for the analysis of large param- 
eter space volumes for coritirmous wave searches. Eaeh step 
only analyzes the regions in paxameter space that have not 
been discarded by any of the previous steps. 



into smaller segments and analyze each segment coher- 
ently. The results of this coherent analysis of the different 
segments are combined incoherently. The output of the 
incoherent step is a set of possible pulsar candidates. If 
necessary, acquire fresh data and repeat the above proce- 
dure analyzing only the candidates selected by the pre- 
vious step. Once this procedure has been iterated the 
desired number of times and the number of candidates 
in parameter space is small enough, the candidates are 
analyzed by using the entire data stream coherently. The 
final output of the search is, of course, either a detection 
or an upper limit. 

The exact number of times the incoherent step must 
be repeated and the thresholds that one must set at each 
stage arc decided by optimizing the sensitivity subject 
to the obvious constraints on the desired signal strength 
we wish to detect, the desired confidence level and the 
amount of total data available. Preliminary investiga- 
tions of this optimization are reported in [9] and more 
detailed results will be presented elsewhere [26]. 

Hierarchical searches like this are typically effective 
only when looking for signals that, in the final coher- 
ent search over the whole data set, have relatively high 
signal-to-noise ratio. The method only works if the in- 
coherent step succeeds in reducing the number of points 
in parameter space that one must search over. A sig- 
nal that is only, say, at two-sigma in the final step will 
be too weak in the initial shorter coherent transforms to 
be selected by any criterion that would eliminate other 
("pure noise") parameter points. Remarkably, this does 
not actually reduce the sensitivity of a hierarchical search 
by much over the hypothetical fully coherent search we 
described above. The reason is that the size of the pa- 



rameter space is so large that, even in a fully coherent 
search, signals must be unusually strong in order to be 
detected with enough significance to be recognized. In 
our case, a fully coherent signal-to- noise ratio of 10 or 
more is needed for a significant detection over a period of 
several months, and we will see in equation (5.35) below 
and the subsequent discussion that our incoherent meth- 
ods do worse than this by factors of between 2 and 5, 
while permitting much larger regions of parameter space 
to be surveyed. 

Finally we mention one important detail, namely, the 
nature of the coherent analysis of each data segment. In 
this paper we consider two possible alternatives. The first 
alternative is just to use the Fourier transform of data 
segments that are so short that no freqTiency modiilation 
or spin-down is measurable. These transforms are called 
SFTs (Short time base-line Fourier Transform) and may 
represent up to 30 minutes of data. The candidates for 
the incoherent step are selected based on the normalized 
SFT power, i.e. on the power divided by the noise floor 
estimate. 

If longer coherent stages are required for better sensi- 
tivity, then one must use demodulated data, i.e. remove 
the effects of Earth's spin and orbital motion and also of 
the pulsar spin-down. This demodulation must be done 
separately for different regions of the sky and spin-down 
parameter space, but it also brings in other parameters, 
such as the polarization angle f/', because of the effects of 
amplitude modulation. These extra parameters, which 
are not part of our Hough-transform search space, can 
be eliminated by requiring the coherent stage to produce 
the J^-statistic described in [7] and used in [6] for ana- 
lyzing the data from the first science runs of the LIGO 
and GEO detectors. In this case, we would select fre- 
quency bins based on the value of the J^-statistic. The 
search based on SFTs will be called the non- demodulated 
search and is described in sections IV and V. The search 
using the JF-statistic is the search with demodulated data 
and is described in section VI. 



III. THE HOUGH TRANSFORM 

As mentioned in the introduction, the Hough trans- 
form is a robust parameter estimator for patterns in dig- 
ital images. It can be used to identify the parameter(s) 
of a curve which best fits a set of given points. In the last 
two decades, the Hough transform has become a standard 
tool in the domain of artificial vision for the recognition 
of patterns that can be parameterized like straight lines, 
polynomials, circles, etc. 

For our purposes, a pattern is a collection C of smooth 
hypersurfaces [32] in some differentiable manifold M. As- 
sume that there is a manifold E of parameters which 
describes elements of C; i.e. there exists a function 
/ : S — > C providing a one-one association between points 
in S and elements of C. 

A simple example is the case when M is ffi.^ with coor- 



5 



dinatcs x and y, and C is the collection of straight lines 
in this {x, y) plane. Since all straight lines arc described 
by an equation of the form y = mx + c (the master equa- 
tion), the parameter space E is also M^, with coordinates 
(m, c) - the slope and the y- intercept of the straight lines. 
The function / maps the point (m, c) to the straight line 
y = mx + c. The relevant example for our purposes is the 
case when the manifold E represents the pulsar param- 
eters ^ = (/(o), {./(n)}j n) and M is the time-frequency 
plane. The pattern in M is described by the Doppler 
shift formula of equation (2.4). Each value of ^ deter- 
mines the frequcnc;y evolution f{t) and thus determines 
a curve in the time-frequency plane. 

Given a set of observations {xi} with each Xi belonging 
to M, we ask if there is an underlying pattern describing 
these points and whether this pattern is described by a 
hypersurface belonging to C. Consider first the idealized 
case when there is no noise and the points {xi} actually 
do follow the pattern and lie on one single hypersurface 
belonging to C corresponding to the parameter value jl € 
S. How would we go about finding /t if we were given 
the collection {xj}? For every Xi, the idea is to first find 
the set of points Ui in parameter space consistcint with 
Xi; the true parameter value p, must certainly lie within 
this set. In the straight line example, all the lines passing 
through the observed point would be consistent with that 
observation. Repeating this for every observation Xi , we 
obtain a collection of subsets {Hi}. The true parameter 
value jj, must lie in each Ui and therefore it must also lie 
in the intersection 

Aefjz^i. (3.1) 

i 

See figure 3. If fc is the dimensionality of S, then we need 

at least k different .x^'s in order to ensure that fi can be 
found uniquely. Thus, in this idealized noiseless case, we 
would need only two observations to detect a straight 
line. Similarly for the pulsar case, equation (2.4) is the 
master equation and if we were searching for s spin-down 
parameters, we would need only 3 + s observations to 
determine the pulsar parameters. This is, of course, not 
true when noise is present. 

In realistic situations, the presence of noise will ensure 
that, in general, there is no point which is consistent 
with all the Xi's, in other words, (liUi is the empty set. 
In this case we proceed as follows: to each /i G S, assign 
an integer n(/x) (the number count), which is equal to the 
number of l/j's which contain fx. The result is then a his- 
togram in parameter space. This procedure, which maps 
a set of observations to a histogram in parameter space, 
will be called the Hough transform. The best candidate 
for the true parameter fi is then the point at which the 
number count is maximal. Alternatively, we could set an 
appropriate threshold riti, on the number count and select 
all points in S at which the number count exceeds ritn- 
These selected parameter space points would be candi- 
dates for a possible detection and, if we were performing 
a multi-stage hierarchical search, would be further ana- 




FIG. 3: A schematic depiction of tlie Hough transform in the 
absence of noise. The top figure shows the parameter space 
S and the space of observations M. The space of expected 
patterns in is a set C of hypersurfaces in M. The function 
/ : E — > C provides a one-one correspondence between E 
and C. The lower figure shows the Hough transform itself: 
Every observation Xi is mapped via the Hough transform into 
a iiypcrsurfacc Ui in parameter space wliicli is consistent with 
the observation. The intersection of all the Ui's contains the 
true source parameter ft. 



lyzed in the next step. 

In real experiments, we cannot perform a parameter 
space search with infinite resolution. Therefore we need 
to consider the discrete case when we have a finite reso- 
lution for the observations and also a grid on parameter 
space. In this case, observations correspond to pixels in 
M. The general procedure is essentially the same as in 
the discrete case and is depicted schematically in figure 4: 
we look for pixels in parameter space which are consistent 
with the observations. There is, however, one technical 
difference namely, since each observation is an extended 
region in M, the points in parameter space consistent 
with this observation do not constitute a sharp hyper- 
surface Ui . Each pixel instead gives a region Ui bounded 
by two such hypersurfaces. Given such a region, we can 
then select pixels in parameter space. Since a pixel in 
parameter space might intersect more than one Ui, we 
need an unambiguous criterion to select pixels in param- 
eter space in order to ensure that each pixel gets selected 
at most once by an observation. Given such a criterion, 
we can continue the earlier strategy and construct a his- 
togram in parameter space by assigning a number count 
to each pixel in parameter space. The pixel with the 
largest number count is our best candidate for a detec- 
tion. 
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FIG. 4: A schematic view of the Hough transform for the dis- 
crete case. An observation consists of a pixel in M which goes 
over to the region enclosed between the dotted lines under the 
Hough transform. This in turn leads to a selection of pixels 
in parameter space. The shaded pixels are the ones which get 
selected and are the ones consistent with the observation. 



IV. THE HOUGH TRANSFORM WITH 
NON-DEMODULATED DATA 

The steps involved in a single incoherent stage of the 
search are outlined in figure 5. In this search, one starts 
by breaking up the input data of duration T<,ba into 
N segments each with a duration of T^oh, which would 
be equal to T„i,jN if there were no gaps in the data. 
Except for precisely two exceptions, namely equations 
(5.35) and (6.41), all the equations in this paper will be 
valid even in the presence of gaps; we shall not assume 
Taah = T^bs/N. This is important because in practice, the 
real data stream will inevitably have gaps in it represent- 
ing times when the detector is not in lock or the data is 
not reliable. 

The next step is to compute the Fourier transform of 
each data segment to obtain N SFTs. Select frequency 
bins in each SFT by setting a threshold on the normalized 
power spectrum. This produces a distribution of points 
in the time-frequency plane — the manifold M — most 
of which are noise but some excess of which are hopefully 
present along one or more signal patterns given by equa- 
tion (2.4). Having selected points in the time-frequency 
plane, go through the Hough transform algorithm to ob- 
tain the Hough map, i.e. the histogram, in parameter 
space S. The details follow. 



A. Notation and conventions 

We assume that the N different data segments have 
the same time duration. Label the different segments by 
a = 0,1... {N — 1) and denote the start time of each 
segment by ta which will often be called the timestamp 
of the a}^ data segment. Let each segment consist of M 
data points. 

Let us now focus on the a*'* data segment which covers 
the time interval \ta,ta + Tcoh]- Let x{t) be the detector 
output which is sampled at times tj = ta + j^t with 
i — 0, 1, . . . {M — 1). Here the data segment has been 
subdivided into M sub-segments with the times tj de- 
fined to be at the start of each sub-segment; this implies 
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FIG. 5: A single stage of a hierarchical continuous wave search 
involving the Hough transform. The starting point is to break 
up the data with total observation time Tobe into segments 
and to compute the Fourier transform of each segment. The 
next step is to select frequency bins from each SFT by setting 
a threshold on the normalized power spectrum and use the 
selected frequency bins to construct a Hough map. The out- 
put is then a set of candidates in parameter space obtained 
by setting a threshold on the Hough number count. 



Ai = T^^^/M . Denote the sequence of data points thus 
obtained by {xj} where Xj = x{tj). 

Our convention for the Discrete Fourier Transform 
(DFT) of {xj} is 

M-l 

ife = At ^ Xj-e-""^''/^'^ (4.1) 

where /fc = 0, 1 . . . (M - 1). For < A; < [M/2J , the 
frequency index k corresponds to a physical frequency of 
/fc — k/T^ah with [.J denoting the integer part of a given 
real number. The values [M/2J < k < M — 1 correspond 
to negative frequencies given by fk^{k — M)/T„„h- 

The detector output x{t) at any time t is the sum of 
noise n(t) and a possible gravitational wave signal h(t) 
of known form: 

x{t) = n{t) + h{t) . (4.2) 

In the remainder of this paper, unless otherwise stated, 
the stochastic process n{t) is assumed to be stationary 
and Gaussian with zero mean. 

In the continuous case, when the observation time is 
infinite, the single-sided power spectral density (PSD) 
Sn{f) for / > is defined as the Fourier transform of the 
auto-correlation function: 

/oo 
{n{t)n{Q))e-^^'f'dt (4.3) 
-OQ 

where (■) denotes the ensemble average. 
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The normalized power is a dimensionless quantity de- 
fined as 



Pk 



(4.4) 



It can be shown that (|nfc| ) is related to the PSD: 

(Kl^) « ^sMk) = ^sMk) . (4.5) 



Thus: 



Pk 



2\xk 



^. |2 



T^ahSnifk) 



(4.6) 



NaturaUy, the PSD must be estimated in a way that is 
not biased by any signal power that may be present. 



B. Implementation 

The implementation choices we present here mostly 
correspond to those that have been implemented in the 
Hough analysis code which is publicly available as part 
of the LIGO Algorithms Library (LAL) [16], and will 
be used to analyze the data from the GEO and LIGO 
detectors. 

Restriction on T„,,; For non-demodulated data, the 
coherent integration time T^oh, i-e. the time-baseline of 
the SFTs, cannot be arbitrarily large. This restriction 
comes about because we woiild like the signal power to 
be concentrated in half a frequency bin but the signal 
frequency is changing in time due to the Doppler modu- 
lation and also due to the spin-down of the star. If / is 
the time-derivative of the signal frequency at any given 
time, in order for the signal not to shift by more than 
half a frequency bin, we must have |/|Tj,oi, < (2Teoh)~^, 
i.e. 



2|/l„ 



(4.7) 



where by |/|„ax we mean the maximum possible value of 
I /I for all allowed values of the shape parameters ^. The 
time variation of f{t) is given by equation (2.4) and is 
due to two effects: the spin-down of the star, and the 
Doppler modulation due to the Earth's motion. We shall 
assume that the Doppler modulation is the dominant ef- 
fect [33] . Thus we can estimate / by keeping / fixed and 
differentiating v(t) in equation (2.4): 



/ 



n < - 

c at c 



dv 
'dt 



(4.8) 



The important contribution to the acceleration dv/dt is 
from the daily rotation of the Earth: 



I j max 



/ 



Re 



c ■ T2 



(4.9) 



where Vg is the magnitude of the velocity of Earth around 
its axis, Te the length of a day and Re the radius of Earth. 
Substituting numerical values we get 



T„h < 50 min x 



/ 500 Hz 



In this paper, we shall mostly use T^, 
canonical reference value. 



(4.10) 



30min as the 



Selecting frequency bins: The simplest method of 
selecting frequency bins is to set a threshold pth on pk] 
i.e. we select the /c*'' frequency bin if pk > pth and reject 
it otherwise. Alternatively [24, 27], we could impose 
additional conditions such as requiring that pk > Pk+i 
and Pk > Pk-i, i-e. the fc*'' bin is selected if pk exceeds 
the threshold and is, in addition, a local maxima. 
This can be extended further by including more than 
just the two neighboring frequency bins. While it is 
relatively easy to investigate these alternate strategies 
for non-demodulated data, the analysis becomes more 
complicated for demodulated data. Furthermore, while 
these alternate methods might be more robust against 
spectral disturbances, the analysis of the statistics 
follows the same general scheme and the results are not 
qualitatively different. Thus, for the purposes of this 
paper, we will describe only the simple thresholding 
scheme for selecting frequency bins. The optimal choice 
of the threshold pth is described below in subsection VB. 

Solving the master equation: As discussed in the pre- 
vious section, to perform the Hough transform, we must 
find all the points in parameter space which are consistent 
with a given observation. In this case, the observation is 
a frequency fk selected using a threshold pti, in say, the 
gth gprj, corresponding to a timestamp ta- This corre- 
sponds to a frequency bin {fk — ^6f, fk + ^Sf) where 
5f = T~l is the frequency resolution of the SFT. The 
parameters ^ of the signal are the frequency, spin-down 
parameters, and the sky-positions: ^ = (/(o), {/(n)}> n). 
Corresponding to {fk — \5f, fk + \5f), we must find all 
the possible values of ^ which satisfy the master equation 
(2.4). 

To understand this better, let us first fix the values of 
the frequency /(g) and the spin-down parameters {/(„)} 

so that f{t) is also fixed. Ignore, for the moment, the 
frequency resolution 5f. From equation (2.4), we see 
that all the values of n consistent with the observation 
f{t) must satisfy 



COS( 



r{t) 



c f{t)-f{t) 



"(t) v{t) f{t) 



(4.11) 



where </> is the angle between v(i) and n. This implies 
that the angle (j) must be constant; in other words, the set 

of sky positions consistent with an observation f{t) form 
a circle in the celestial sphere centered on the vector v 



FIG. 6: The set of sky positions consistent with a given fre- 
quency bin at a given time correspond to annuli on the celes- 
tial sphere. These annuli are centered on the velocity vector v, 
they are thin when perpendicular to v and thick when nearly 
parallel. The circle with with the label n = corresponds to 
/ = /■ 



(see figure 6) [34]. If the frequency f{t) is smeared over 
a frequency bin (/^ — ^Sf,fk + |^/), the set of points 
consistent with an observation must correspond to an 
annulus the width 6ip of which is easily calculated using 
equation (4.11): 



5f 



V /sin( 



(4.12) 



The annuli arc very thick at points where sin(j) is small, 
i.e. when n is almost parallel or anti-parallel to v{t) and 
very thin when perpendicular. This is depicted schemati- 
cally in figure 6. The circles on the celestial sphere are la- 
beled by an integer n such that the frequency / = f+nSf 
corresponds to the angle ^„ given by 



COS( 



ncdf 



(4.13) 



The lower limit on the width of the annuli is provided 
by setting = 7r/2 in equation (4.12): 



= 4.8x 10"^radx 



lhr\ /500Hz\ /lO"'' 



(4.14) 



^coh /V / 



v/c 



The upper limit on the annuli width ((5(^)„„ is found by 
setting sin(/> w w (5?!>)max which gives 




.r.-2 J / Ihr V / SOOHz V / 10"^ 

= 7.3 X 10^ rad x — — — — 

\T.oJ \ f J \v/c 

Therefore, the thick annuli arc about 10 times thicker 
than the thin ones. Different frequency bins selected at 



FIG. 7: Two intersecting annuli. The two timestamps to and 
tb are sufficiently different from each other so that the ve- 
locities v(ta) and v{tb) arc not parallel to each other. This 
causes the annuli constructed at different timestamps to in- 
tersect. The shaded region is the intersection and there is a 
corresponding region (not shown) on the far side of the sphere. 



the same time will correspond to non-intersecting annuli 

as shown in figure 6. However, for frequency bins selected 
from SFTs at different time stamps, say ta and t^, the 
annuli will usually intersect because the velocity vectors 
v(ia) and v(tf,) will not, in general, be parallel to each 
other; see figure 7. p 

Resolution in the space of sky-positions: In order 
to search for pulsar signals in a given portion of the sky, 
we must choose a tiling for the sky patch. Given the cal- 
culation of the annuli width above, we choose the pixel 
size 66 of the grid to be some fraction, say at most half, of 
the width ((5(^)„i„ of the thinnest annulus. While this ed- 
ucated guess for the pixel size is sufficient for the purposes 
of this paper, the correct choice of pixel size in the sky 
patch, and also in the entire parameter space, should use 
the parameter space metric introduced in [28] . The anal- 
ysis of this metric for the Hough search will be presented 
elsewhere, and for now we shall simply use 69 = ^ (^?5')min- 
Having selected an annulus and having chosen a tiling 
on our sky-patch, we now need a criterion for selecting 
a pixel if it intersects an annulus. Our criterion is to 
select a pixel if its center lies within an annulus. Under 
such a criterion, a given pixel can then be selected by 
at most one annulus and the pixels selected by all the 
annuli together will exactly cover the sphere. 

Resolution in the space of spin-down parameters: 

In the absence of a proper analysis of the parameter space 
metric, we shall just use the obvious estimate for the 
resolution ^/(„): 

SfM=n\^. (4.16) 

As an example, for the first spin-down parameter: 

,^ in„ ,^ 30days 1800s 
= (2.1 X 10-^°Hz/s) X — . (4.17) 



'- obs 



T 
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We now need to choose the range of vahies —f^^^ < 
f(n) < f(n) ^"^"^ largest number of spin-down param- 
eters s„ax to be searched over. Assuming that the pul- 
sar's frequency evolution is well-represented by a Taylor 
expansion, we get 



(4.18) 



where r is the age of the pulsar and /^^x is the largest in- 
trinsic frequency that we search over. We include the n*'* 
spin-down parameter in our search only if the resolution 
defined by equation (4.16) is not too coarse compared to 



/max , 
(n) • 



(4.19) 



i5/(n) < /(Tf • 

Since T^t^ <c t, f^^^ decreases with increasing n much 
faster than 5f(ny, this implies that there must exist a 
value s„,ax such that equation (4.19) is satisfied for all 
n < s„ax and is violated for all n > s„,„. Any spin-down 
parameter of order greater than ,s„„„ does not signifi- 
cantly affect the result of the Hough transform. As an 
example, if we wish to search for pulsars whose age is at 
least T = 40yrs, then for /„,,„ = lOOOHz, it is easy to 
check that we get s„ax = 3. In other words, to look for 
pulsars which are as young as 40yr, we must include at 
least 3 spin-down parameters in our search. 

On the other hand, in some cases, computational re- 
quirements might dictate that we can only search over, 
say, one spin-down parameter. This automatically sets a 
lower limit on the age of the pulsar that we can search 
over because then the second spin-down parameter must 
satisiy (5/(2) > /^r which leads to 



T > 155yr x 



1/2 



SOdays V lOOOHz 1800s 



(4.20) 



Finally, the finite length of T^oh itself leads to a lower 

bound on r. If is too large, then the signal power 
may no longer be concentrated in a single frequency bin 
and the assumption of neglecting spin-down parameters 
which was used to derive equation (4.10) will no longer 
be valid. To be certain that the spin-down will not cause 
the signal to move by more than half a frequency bin, we 
must have f^^T^i, < n\Sf/2 which implies 



r > 



o f q-'n+l 

max-*- coh 



1/n 



n! 



The most stringent limit is obtained for n = 1: 



T > 103yr X 



tax / T^.. 



lOOOHz V 1800s 



(4.21) 



(4.22) 



This restriction will not be present if we use demodulated 
data as input for the Hough transform. 
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FIG. 8: A Partial Hough Map (PHM) is a histogram in the 
(a, 5) plane constructed from all the frequencies selected at 
a given time and for a given value of the instantaneous fre- 
quency fo- A total Hough map is obtained by summing over 
the appropriate Hough maps. The PHMs to be summed over 
are determined by the choice of spin-down parameters which 
give a trajectory in the time-frequency plane. For example, a 
single spin down parameter will give a straight line as shown 
is the figure while two spin-down parameters will lead to a 
parabola. 



Partial and total Hough maps: As described above, 
for a given frequency bin selected at a given time-stamp 
and for a given value of the instantaneous frequency /, 
we can find the set of sky locations which arc consis- 
tent with the master equation (2.4). In other words, 
every pixel in the sky-patch either gets selected or re- 
jected and this gives a histogram in the (a, S) plane con- 
sisting of ones or zeros; a and 6 are coordinates on the 
sky-patch. Such a collection of ones and zeros on the 
sky-patch is called a Partial Hough Map (PHM). The 
number of PHMs required at any given time depends on 
the frequency band A/f, that one is searching over and is 
given by Ah/Sf = T^^Afb- 

Given a set of PHM's for every time interval, and 
given a set of spin-down parameters that one wishes to 
search for, the Total Hough Map (THM) is obtained by 
summing the appropriate partial Hough maps. To see 
how this comes about, consider the case when we are 
searching for some spin-down parameters {/(„)} with 
n = 1, 2, • • • . The instantaneous frequency changes with 
time according to equation (2.5); ignore the Ar term in 
this equation [35]. This can be viewed as a trajectory in 
the time-frequency plane. A single spin-down parameter 
will give a straight line, two spin-down parameters a 
parabola and so on. Thus for each time-stamp ta, we 
can find the appropriate PHM by looking at which 
frequency bin this trajectory intersects (see figure 8). 
For a given choice of spin-down parameters, the THM 
is obtained by summing over the appropriate PHMs. 
Repeating this for every set of frequency and spin-down 
parameters we wish to search over, we obtain a number 
of THMs and the collection of all these THMs represent 
our final histogram in parameter space. 
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Look up tables: The procedure described thus far is, in 
principle, enough to produce a complete Hough map in 
parameter space. However, it is possible to enormously 
reduce the computational cost by using Look Up Tahles 
(LUTs) which we now describe. Assume that we have 
managed to find all the annuli for a given time-stamp 
ta and for a given search frequency /. To construct the 
PHM for ta and /, we just need to select the appropriate 
annuli out of all the ones that we have found. Very im- 
portantly, it turns out that in most cases, the annuli are 
relatively insensitive to changes in / and can therefore 
be re- used a large number of times. 

To see this, look at how the solutions of equation (2.4) 
depend on the search frequency /. Wc want to calculate 
the maximum number n of frequency bins that / can be 
changed by so that the annuli change by only a fraction 
r of the quantity ((5(/))inin defined in equation (4.15). As 
discussed earlier, if we restrict ourselves to discrete fre- 
quencies, the annuli corresponding to every given value 
of / are parameterized by an integer n according to equa- 
tion (4.13). For a fixed value of n, by how much does 0„ 
change when / is varied? To answer this, differentiate 
equation (4.13) with respect to /: 



let Np be the number of templates in the space of sky- 
locations and spin-down parameters; T^^^^Afh is the num- 
ber of frequency bins being searched over. 

A naive implementation of the Hough transform will 
require, for every point in parameter space, to iden- 
tify first the corresponding pattern in the time- frequency 
plane and then add N integers (zeros or ones) to obtain 
the final number count, where N is the number of data 
segments. If we use LUTs, which are valid for a large 
number of search frequencies, their computational cost 
in the search becomes negligible, i.e. the cost of find- 
ing the patterns is negligible, and the number of fioat- 
ing point operations required is thus Co « T^^i^AfbNpN. 
This calculation can be organized much more efficiently 
if we perform a search on many sky locations at once. In 
this case, if we know the locations of all the annuli, for 
every chosen frequency bin, we mark the corresponding 
annulus and in the end, combine all the annuli thus se- 
lected to get the final number count. The exact savings 
in computational cost due to this strategy are implemen- 
tation dependent, but are typically better by a factor of 

5 when compared to the value Cq mentioned above. 
This factor is related to the number of frequency bins 
selected from each SFT. 



-J- = ——-SirKpn d(j)n = f d(j)n ■ 



(4.23) 



STATISTICAL PROPERTIES OF THE 
HOUGH MAPS 



Set d(t)n = r{6(p)^i„ and df = nSf = k/T^oh to obtain 



rc 

K = — tan ( 

V 



V \ in? 



(4.24) 



where no = vf/{cSf). Consider separately the two 
regimes when <pn ~ 7r/2 (i.e. n ~ 0) and <pn ~ 0,7r 
(i.e. n ~ ±rio). When 0„ = 7r/2, then k is infinite which 
indicates that a LUT is excellent in this regime. On the 
other hand, k = for = 0,7r. However, since the 
resolution in 4> is finite, instead of 0„ = 0, it is more ap- 
propriate to take the worst case scenario as 0„ = (50)„i„ 
so that 



/ 500 Hz 

— (66) . = 40r — 

V V / 



(4.25) 



Thus, in this worst case scenario, for a frequency of 500 
Hz and a tolerance of r = 0.1, the LUT will be valid for 
4 frequency bins. Furthermore, due to the presence of 
the function tant^ in equation (4.24), n increases rapidly 
with increasing (i.e. decreasing n). As an example, 
take T,oh = 1800s, / = 500 Hz, and v/c = IQ'^ so that 
no = 90. Then, even for n = 89, we get k = 1500r; 
thus with say r = 0.1, the LUT is valid for about 150 
frequency bins. 

The main point of using the look up tables and partial 
Hough maps is to reduce the computational costs. As- 
sume that we are searching over a frequency band A/(,, 



This section is divided into three parts: The prob- 
ability distribution of the number counts is calculated 
in subsection VA, subsection VB optimizes the various 
thresholds and subsection VC estimates the sensitivity 
of the Hough search. 



A. The number count distribution 

The frequency bins that are fed into the Hough trans- 
form are the ones such that their normalized power pk 
defined in equation (4.4) exceeds a threshold pth- As- 
suming that the noise is stationary, has zero mean, and 
is Gaussian, from equation (4.4), we get 



2pfe = Z-i^+Z2 



where 



Zl 



V2Ri 



e\Xk\ 



and 



Z2 



m\Xk\ 



(5.1) 



(5.2) 



As before, the detector output aifc is the sum of noise and 
a possible signal: Xk = fik + hk- Assuming that Re[nfc] 
and Im[nfe] are independent random variables with equal 
variance, it is easy to show that their variance must be 
equal to (|nfcp)/2. Therefore, taking the noise to be 
Gaussian, it follows that the random variables zi and 
Z2 are normally distributed and have unit variance (but 
non-zero mean). Thus 2pk must be distributed according 
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to a non-central distribution with 2 degrees of freedom 
with non-centrahty parameter A^: 

A, = (E[.,])^ + (Eh])^ = . (5.3) 

Thus the distribution of pk is 

p{pk\Xk) = 2x\2pk\2,Xk) 

= exp (^-pk - Ia{\/2XkPk) (5.4) 

where Iq is the modified Bessel's function of zeroth or- 
der. As expected, p{pk\Xk) reduces to the exponential 
distribution in the absence of a signal (when A = 0). 

The mean and variance for this distribution are respec- 
tively 

Ek] = l + Y and a2[pfe]=l + Afc. (5.5) 

The probability t] that a given frequency bin is selected 
is 

V{P..\X)= / p{p\X)dp (5.6) 

where we have dropped the subscript k for notational 
simplicity; it is understood that p and A always refer to 
one of the Fourier frequency bins. The false alarm and 
false dismissal probabilities for the frequency bin selec- 
tion are respectively 

/■OO 

a(pth) = / pip\0)dp = e-P'\ (5.7) 

/3(Ah|A) = l-v{p,^\X)= r\{p\X)dp. (5.8) 

Jo 

Clearly, r] = a when no signal is present and r] becomes 
larger when the signal becomes stronger and 77 1 when 
A ^ 00. Figure 9 shows r]{pti,\X) as a function of the non- 
centrality parameter A for two different values of pth- For 
small A: 

rj{p,,\X) = a{l + ^X + 0{X')} . (5.9) 

This expansion will be very useful when we restrict our- 
selves to the case of small signals. 

In the presence of a signal, the non-centrality parame- 
ter Afc is not constant across different SFTs. The reason 
for this is two- fold: First, the noise may be significantly 
non-stationary. Secondly, and more fundamentally, the 
observed signal power changes because of the ampli- 
tude modulation of the signal caused by the non-uniform 
antenna pattern of the detector. Therefore, the detection 
probability r? changes across SFTs. In what follows, we 
shall neglect this effect and take A and r] to be constant 
for different SFTs. 
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FIG. 9: Plot the detection probability ?7(pth|A) as a function 
of A for pth = 1.6 and 2.6. 

Under this assumption, the probability of measuring a 
number count n in a pixel of a Hough map constructed 
from N SFTs is given by the binomial distribution: 

p(n|p,„A)=(^)r?"(l-r?)^-". (5.10) 

The mean and variance of the number count are respec- 
tively 

n = Nr} and cr^ = Nr}{l - rf) . (5.11) 
In the absence of a signal, rj = a so that 

p(nK,0)= (^)a"(l-a)^-". (5.12) 

Candidates for detection or for further analysis are se- 
lected by setting a threshold riji, on the number count. 
Based on this, we can define the false alarm and false 
dismissal rates respectively as: 

N 

aH{n,i„p,h,N) = ^ p(n|pth,0), (5.13) 

n=nth 
"th-l 

X,N) = J2 f("lA'>,A). (5.14) 

n=0 

These two quantities determine the significance and the 
sensitivity of the Hough search and will play an impor- 
tant role in the rest of this paper. 

B. Optimal choice of the thresholds 

In order to carry out the Hough search, we have to 
set two thresholds: the threshold p^^ on the normalized 
power and the threshold rith on the number count. 

The value of rith is determined by the false alarm rate 
that depends on the number of candidates that we 
can feasibly follow up. 
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The value of pth is chosen in such as way so as to 
make the search as powerful as possible. We present two 
criteria that yield the same result for small signals and 
for large N. 

Maximizing the critical ratio: For the Hough num- 
ber count, we can define a random variable called the 
critical ratio as follows 



* = 



n — Na 
^/Na{l - a) 



(5.15) 



This quantity is a measure of the "significance" of a mea- 
sured value n with respect to the expected value Na in 
the absence of any signal, weighted by the expected fluc- 
tuations of the noise. In the presence of a signal, the 
expected value of the critical ratio is 



*(r?,a) 



Nt] - Na 



(5.16) 



Recall that rj and a depend on the threshold p^^ according 
to equations (5.7) and (5.8) respectively. Thus, our cri- 
terion for choosing the threshold is to maximize ^{r],a) 
with respect to pti,- In the case of small signals where 
r] ~ a(l -|- pthA/2), the condition 



9* 
dpth 



= 



(5.17) 



leads to 



lna = 2(a-l) (5.18) 
which yields ~ 1-6 or equivalently, a ~ 0.20. 

The Neyman-Pearson criterion: An alternative 

method of choosing pti, is based on the Neyman-Pearson 
criterion which tells us to minimize the false dismissal 
rate Ph for a given value of the false alarm rate. 
For weak signals, this requirement uniquely determines 
Pth and, as we shall see, this agrees with the previous 
criterion. 

In practice, taking N and A to be fixed parameters, 
this is the procedure: 

i. First choose a value for the largest false alarm rate 

an that we can allow. 

ii. Invert the equation aH{Pt^-iN,nt_i,) < to obtain 

iii. Substitute the value of t;,i, thus obtained in the ex- 

pression for the false dismissal /3i/(nth, Pth, A, iV). 
This gives as a function of (pth, A, N, a^). 

iv. Minimize (3h as a function of Pth- Let p*^^ be the value 

that minimizes (3h- 

V. Using nth(pth, N, a^) derived in step (ii) above, obtain 




FIG. 10: Graph of nth — A'^Q versus the false alarm probability 
a = e""'" for an = 0.01 and N = 2000. The dashed line 
shows the analytic approximation given by equation (5.21). 



This procedure is also applicable if we choose a differ- 
ent method of selecting frequency bins other than simple 
thresholding, such as, for example the peak selection cri- 
terion mentioned towards the end of subsection IV A. 

The results of the optimization procedure described 
above are shown in figures 10, 11 and 12. Figure 10 
shows the value of the number count threshold n^^ ob- 
tained as described in step (ii). In this figure, instead 
of pth, we have chosen the false alarm rate a = e"''"' as 
the independent variable; a is the false alarm rate for 
selecting frequency bins and is not to be confused with 
an- Figure 10 also shows an analytic approximation to 
rith obtained below in equation (5.21). Using this result 
for nth, figure 11 shows as a function of a = e"''*''. 
The optimal choice p*^^ of Pth is when I3h is a minimum 
and, for small signals, this happens at p*^ s:^ 1.6 which 
corresponds to a* := e"''"- ~ 0.20 . Finally, figure 12 
shows the minimum value of obtained by this opti- 
mization as a function of the signal strength A and for 
two different values of N. 

The Gaussian approximation: To better understand 

the statistics, it is useful to carry out the above steps an- 
alytically by taking n to be a continuous variable and by 
approximating the binomial distribution by a Gaussian 
with the appropriate mean and variance: 



p(n|pth,A) 



\/27ra2 



(5.19) 



This is a very good approximation when is large and 
r\ is not too close to or 1. If n is chosen to lie within 
[0,A], the distribution is properly normalized only ap- 
proximately. For simplicity, in what follows we shall take 
the range of n to be (— oo, -l-oo); this is an excellent ap- 
proximation if the above assumptions on N and 77 hold. 

With the approximations given above, we can rewrite 
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FIG. 11: The first figure shows the Hough false dismissal rate 

as a function of the false alarm rate a — e"''"' for a non- 
centrality parameter A — 0.10 (upper curve) and A = 0.20 
(lower curve). Both curves correspond to oh = 0.01 and 
A'^ = 1000. The minimum values of Ph for the two curves are 
approximately 0.84 and 0.41 respectively. Both minima occur 
at Q — 0.20 approximately. This corresponds to a threshold 
of pth = 1.6 on the normalized power statistic. The bottom 
figure shows the approximation to /3h using equations (5.22) 
and (5.21) with the same paxameters as in the first figure. 



the equation an = as 

p{n\pt^,0)dn = ■ 



I 



(5.20) 



The solution to this equation can be rewritten in terms 
of the complementary error function: 



n,H(Ah, N, a*H) = Na+ ^2Na{l - a) erfc-^(2a^) . 

. ^^-^^^ 

As shown in figure 10, this is a very good approximation 
to the actual value of rith obtained from the binomial 
distribution. 

The expression for I3h is similarly rewritten as 



I3h = ;^erfc 



/ Nrj - rtth 
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FIG. 12: Minimum value of Ph as a function of the non- 
centrality parameter A for oh = 0.01 and for N = 1000 and 
2000. As expected, a larger value of N typically leads to a 
smaller value of (Sh- 



As figure 11 shows, this too is a very good approximation 
to Ph obtained using the binomial distribution. 
In step (iv), we find p*^ such that 



dfiH 



dpa 







(5.22) 



In the case of small signals where rj « a{l + pthA/2), the 
solution to the above equation becomes independent of 
A, and it is also independent of N when N is large. In 
these limits, the solution is given by 



d 



ePth - 1 



(5.23) 



The solution to this equation is p*^^ « 1.6 and a* = 
e~Pth Ki 0.20. Notice that this equation is equivalent 
to equation (5.17) and furthermore, the functions being 
extremized are rather flat near the extremum. Thus, the 
threshold could be chosen somewhat diff'erently without 
significantly impacting the sensitivity. In particular, the 
threshold can be increased so that fewer frequency bins 
are selected. Depending on the details of the implemen- 
tation, this could lead to a lower computational cost; in 
the framework of a hierarchical search, this will improve 
the overall sensitivity. 

Finally, with the optimal threshold p*^^ at hand, the 
optimal threshold n*^ on the number count is obtained 



by substituting p^, = p*^ in equation (5.21): 
<h = n,4pl,N,a*H) 



(5.24) 



= Na* + ^/mct*{^^^c^ eTk~\2a*H) . 

This is an important equation because it tells us the num- 
ber count threshold that must be set in order to have a 
given number of follow-up candidates. 
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C. Sensitivity 

In this subsection, we estimate the sensitivity of the 
Hough search, i.e. we answer the foUowing question: for 
given values and of the false alarm an and false 
dismissal /3h respectively, what is the smallest value of 
the gravitational wave amplitude ho (see equation (2.9)) 
that would cross the thresholds p^t, and rith? Equivalcntly, 
for a given false alarm rate a^, what is the smallest ho 
which will give a false dismissal rate of at least We 
use the signal model of equation (2.7) and we present 
our final result for the values Uh = Oi*H — ^-Ol ^^"^ 
(3h = l^H — 0-10. The value of 0.01 is meant mainly for il- 
lustration purposes and does not change the results qual- 
itatively. Furthermore, for comparison, equation (2.2) in 
[6] assumes a false alarm of 0.01 and this choice of a*^ 
enables an easier comparison with that result. As far 
as possible, we explicitly retain the factors of in our 
equations substituting numerical values only when nec- 
essary. 

We must first solve the equation 



/3h«,p*„A,7V)=/3^ 



(5.25) 



and obtain A as a function of N; this will yield the desired 
value of /lo- 
in order to simplify the discussion, we shall again ap- 
proximate the binomial distribution by a normal distri- 
bution whose mean n, and variance cr, are respectively 
given by equation (5.11). The false dismissal rate is 



P. 



1 



H 



rfc 



V'2iVr?*(l - rj*) 



(5.26) 



where n*^^ is as given in equation (5.21) and 77* = ?7(/5*i,l'^)- 
Since we are interested in the case of small signals, let 
us approximate rj by only keeping terms of the order of 
A in equation (5.9). Ignoring terms of C(A^), equation 
(5.26) leads to the approximation 

Hh = ^erfc (^-erfc-^(2a^) + ^OaV^.A^ (5.27) 



where 



e 



N 



2a*(l - a*) 



1 - 2a^\ erfc-^(2a^) 



1-a* 



2a* 



(5.28) 
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FIG. 13: Graphs of (3h as a function of A for different values 
of N and for the three different approximations used. In the 
first panel N = 20, the second panel has N = 2000 and in the 
third panel, A'^ = 20, 000. All graphs arc plotted assuming 
the optimal values for pth and rith- The linear approximation 
is clearly unacceptable for A'^ ^ 10^ but becomes reasonable 
when N ~ 10^ or 10^ and is excellent for N r^W^. The Gaus- 
sian approximation is clearly much better and is good even 
for N = 20. Finally, note that the approximations always 
underestimate the value of (3h- 



Let us summarize our approximation scheme for (3h- The 
first approximation is to take the number count distri- 
bution to be binomial. The second approximation is in 
equation (5.26) which replaces the binomial by a Gaus- 
sian distribution with the appropriate mean and vari- 
ance. The final approximation is in equation (5.27) where 
we have taken A to be small and used a Taylor series in 
powers of A retaining only the linear term. To get a 
feeling for the validity of these approximations, figure 13 
shows graphs of Ph as a function of A for different values 



of N. As the graphs show, we can trust the approxima- 
tions when N ~ 10'^. For smaller values of N, while the 
Gaussian approximation is still reasonable, the linear ap- 
proximation greatly under-estimatcs (Sh for a given value 
of A, i.e. it make the Hough search appear more sensitive 
than it actually is. 

Working with the linear approximation of equation 
(5.27), assuming N to be very large and erfc~^(2aj^) <C 
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N, set I3h = Ph ^^'^ solve for A: 



S_ 



8(1 - a*) 9.02 



Na* 



where 



S := evfc-\2a*H) + erfc-^(2/3^) 



(5.29) 



(5.30) 



and to obtain numerical values, we have chosen = 
0.01 and = 0.10. Using the properties of the com- 
plementary error function, it is easy to show that 5 = 
implies that the statistical significance s := 1 — — f}'^ 
also vanishes. Therefore, the quantity <S can be taken to 
be a measure of the statistical significance of the search. 
The value of A obtained in equation (5.29) gives us the 
strength of the smallest signal that can be detected by 
the Hough search with a false alarm rate of 1% and a 
false dismissal rate of 10%. 

A graph of A as a function of for small values of 
N is shown in figure (14); this figure shows the results 
using both the linear approximation and the more accu- 
rate binomial distribution. The small N limit requires 
a brief explanation. For small N . the discrete nature of 
n becomes important. In particular, the false alarm an 
defined in equation (5.13) can take only a discrete num- 
ber of values, the smallest of which is (at = N). 
Thus for = 1, it is not possible to reach the desired 1% 
false alarm rate and the best we can do, with 

Pth — 1.0, is 

an = 0.2. To find the value of A which yields /3h = 0.1, 
note that for iV = 1, /3// = 1 - rj. Thus r]= 1-0.1 = 0.9 
which implies A w 8.08; this is the sensitivity of the 
search for A = 1. It corresponds to a false alarm rate 
of 20% and a false dismissal rate of 10%. Similar calcu- 
lations show that the sensitivity becomes worse as A^ is 
increased from 1 to 4 as the corresponding false alarm 
rates become better. It is only at A' = 5 that we can 
choose n^i, < A from which point onwards the sensitivity 
begins to improve. This explains the small A^ behavior 
of figure 14. Similarly, the other discrete jumps in figure 
14 are due to the discrete nature of ajj and requirement 
of keeping it below the 1% level. 

To recast the expression for A directly in terms of the 
signal amplitude, start with equation (5.3): A depends on 
the various pulsar parameters. The relevant quantity for 
the purposes of this subsection is the average of A over 
these parameters. It is quite straightforward to estimate 
this average. First, recall the expression for h{t): 

h(t) = F+(t)A+ cos<I>(t) + Fx(t)^x sin$(t) . (5.31) 

Since T^oh is much lesser than a day (see equation (4.10)) 
we can take F+x to be roughly constant. Similarly, as- 
suming that r^oh is small enoiigh so that the pulsar signal 
does not shift by more than half a frequency bin, we can 
take the signal frequency f{t) to be roughly constant. 
With these approximations, we get 



^ 7r(/-/fe)Jeoh 
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FIG. 14: Graph of the smallest detectable A with the optimal 
thresholds. The dashed curve uses the linear approximation 
of equation (5.27) while the solid curve uses the binomial dis- 
tribution. See text for additional discussion. 



where / is the instantaneous frequency of the signal and 
fk is the central Fourier frequency of the frequency bin 
containing /; / is allowed to lie in the range {fk — 
Sf/2, fk + 5f /2). Now take the square of hk and average 
over time to get the average non-ccntrality parameter for 
all the SFTs and note that the time averages and 
are both 1/5, and the time average of F+F^ vanishes. 
Thus: 



A 



105„ 



(Al+Al) 



sin[7r(/ - fk)T,^h] 

7r(/ - /fe)Teoh 



(5.33) 



Take the amplitudes to be of the form given in equations 
(2.7) and average over cos l G (—1, 1) and over the values 
of the signal frequency f G {fk - Sf/2, fk + 5f/2): 



4 /igr,„, 

25 Sn 



sin^ (ttx) 
{ttx)^ 



dx 



0.7737 X 



4 hlT^^ 

25 Sn 



(5.34) 



We get the following value for the smallest signal that 
can be detected by the Hough search: 



8.54 

Ai/4 



^ = 8.54Ari/4j|^. 



(5.35) 



Here the second equality assumes T^i,^ = NT^^i, which be 
true only if the N different data segments were contigu- 
ous; this is done only to compare this result with equation 
(5.36) below. 

Equation (5.35) is the result we were looking for. This 
tells us that if we wish to detect a signal with a false alarm 
rate of 1% and a false dismissal rate of 10%, the weak- 
est signal that will cross the optimal thresholds is the Hq 
, , given above. The important feature to note is that is 
^ ' ' proportional to N^^'^j^jT^^ while for a coherent search 
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over the whole observation time, the sensitivity is pro- 
portional to l/^/r^bB- In particular, for the same values 
of the false alarm and false dismissal rates as above, the 
sensitivity of a full coherent search directed at around a 
single point in parameter space is given by (see [6]): 



ho = 11.4W-^ 



(5.36) 



This illustrates the loss in sensitivity introduced by com- 
bining the different SFTs incoherently but, of course, 
this is compensated by the lesser computational require- 
ments for the incoherent method. Furthermore, for say 

= 2000, the sensitivity of the Hough search is only 
about a factor of 4.5 worse than a full directed coherent 
search. 

This result helps one to make trade-offs of coherent 
against hierarchical searches. For example, if one is 
searching for a population of objects that is uniformly 
distributed in a plane, such as a population of young pul- 
sars in the Galaxy, then a coherent search of any region 
of parameter space would go 4.5 times deeper than the 
incoherent method with N = 2000. The volume of space 
surveyed in the plane would be 4.5^ = 20 times larger. 
However, if the incoherent method's speed of execution 
allowed it to survey more than 20 times as much parame- 
ter space (including sky area and spin-down range) then 
one would choose the incoherent method. This is indeed 
the case for pulsar searches. 

Finally, equation (5.35) also allows us to estimate the 
astrophysical range of the search. Combining (2.9) and 
(5.35), we get: 



167r2GiVi/4/^^e/2 / t^^^ 



8.54c4 



Sn {fr 



(5.37) 



= 5.8kpcx 



N 



f fr 



ViQ-e; Visoosy 



17000 y VlO^®kg-mV V 500Hz 
lO-^^Hz"^^ ' 



Here the reference values for Teoh and N have been chosen 
such that NT^^-t, k, lyr, and we have taken the detector 
sensitivity to be 10~^^Hz~^^^ at a frequency of 500Hz, 
which is appropriate for the proposed advanced LIGO 
detector [29] . A full directed coherent search does better 
than this by a factor of about ~ (17000)^/'' w 12. For 
the initial LIGO detector, assuming it is 10 times less 
sensitive than advanced LIGO, we see that it does worse 
than this by about a factor of ^ 3. 



VI. HOUGH TRANSFORM WITH 
DEMODULATED DATA 

As equation (4.10) shows, using the Hough transform 

with SFTs as input, necessarily limits the coherent time- 
baseline T^oh, and therefore also the sensitivity of the 



search. To get around this limitation, we need to de- 
modulate each coherent data segment to remove the fre- 
quency drifts caused by the Doppler modulation and the 
spin-down; the only limitation on T^oh is then due to the 
available computational resources. The demodulation 
procedure we use is based on the jF-statistic introduced 
in [7]. The search pipeline is very similar to the pipeline 
shown in figure 5, the only difference being that instead 
of computing the power spectrum, we calculate the T- 
statistic. Subsection VIA provides a brief description of 
the .?^-statistic, the master equation is derived in subsec- 
tion VI B, subsection VIC provides the implementation 
details and subsection VI D describes the statistics. 



A. The .^-'-statistic 

Let x{t) be the calibrated detector output and let h{t) 
the waveform that we are searching for. In order to ex- 
tract the signal h{t) from the noise, the optimal search 
statistic is the likelihood function A defined by 



\nK={x\h) - -{h\h) 
where the inner product (-j-) is defined as 



{x\y) := 2 



i{m{f)+S:%I)y{f) 

Snif) 



d/. 



(6.1) 



(6.2) 



Here, as before, x{f) is the Fourier transform of x{t) and 
Snif) is the one-sided power spectral density. The ex- 
pected waveform h{t) is given by equations (2.6), (2.2) 
and (2.3). The quantity In A is essentially the matched 
filter and is precisely what we should use in order to best 
detect the waveform h{t). However, apart from the pa- 
rameters ^ — (/(o), {/(ti)}j ^i)i 111 A also depends upon the 
other parameters such as the orientation of the pulsar, 
the polarization angle of the wave etc. The .F-statistic 
eliminates these additional variables and enables us to 
search over only the shape parameters ^. 

Following the notation of [7], the dependence of the 
antenna patterns -F+,x on the polarization angle tp are 
given by 

F+{t) = sinC[a(t)cos2V'-|-6(t)sin2V'] (6.3) 
Fx(t) = smC[b{t)cos2tp - a{t)sm2tlj] (6.4) 

where the functions a{t) and b(t) are independent of tp 
and C is the angle between the arms of the detector. If 
we write the phase of the gravitational wave as 



$(t) = 00 + Ht) , 



(6.5) 



then we can always decompose the total waveform h{t) 
in terms of four quadratures as 



h{t) = ^A,hi{t) 



(6.6) 
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where the four amphtudes A, are time independent and 
the hj arc as follows: 



hi{t) 



a{t) cos <f){t) , 
a{t) sm(j){t) , 



h2{t) = b{t) C0S(/)(t) , 
K{t) = h{t)sm(j}{t) . (6.7) 



What this decomposition achieves is a separation of the 
shape parameters ^ from the other pulsar parameters. 
The only unknown parameters in the quadratures hi are 
the shape parameters ^ while the amplitudes Ai are in- 
dependent of ^. The log likelihood function depends 
quadratically on the four amplitudes and we can analyt- 
ically find the maximum likelihood (ML) estimators Ai 
of the amplitudes A^ by solving the set of four coupled 
linear equations 



ainA 



dA, 



= 0, i = l,...,4. 



(6.8) 



The J^-statistic is then defined as the log likelihood ratio 
with the values of the amplitudes A^ replaced by their 
ML estimators: 

The only unknown parameters in the optimal search 
statistic T are the shape parameters ^. 



B. The master equation 

Equation (2.4) describes the expected time-frequency 
pattern when the search statistic is the Fourier trans- 
form; in other words, if the detector output x{t) con- 
tains a true signal with instantaneous frequency f{t), 
then equation (2.4) tells us the value of the observed 
frequency f{t) which would maximize |a;(/)P in the ab- 
sence of noise. If we now use the jF-statistic instead of the 
Fourier transform, the expected time-frequency pattern 
is, as described below, different. 

Before proceeding further, it is useful to distinguish 
the instantaneous frequency /(o) from the other shape 

parameters which we denote by A: ^ = (/(o),A) = 
(/(Q), {/(„•)}, n). Let us the assume that the JF-statistic 

has been computed using the parameters A^ but that 
the detector output consists of a signal with parameters 
(/(o) : ^) ; let us denote the mismatch in the parameters 
by AA := A - A^ and A/ := / - /(o). 

Since the .?^-statistic is maximized when the source pa- 
rameters are equal to the demodulation parameters, it is 
clear that if AA = 0, then the expected time-frequency 
pattern is just a constant frequency / = /(o) • More gen- 
erally, due to the correlations in parameter space, the 
mismatch AA may produce a residual shift in the fre- 
quency A/. This frequency shift is determined by 



dJ^{f, Arf;/(o), A) 
df 



= 0. 



(6.10) 



Expand T in powers of AA and A/ around the point 
(/(o) ) ^) up to second order (repeated indices are summed 
over): 

^(/,Ad;/(o),A) = .F(/(o),A) + Aoo(A/)2 (6.11) 
+ ^oiAAiA/-F AyAAjAAj . 

The linear terms do not appear because is maximized 
when A/ = and AA = 0. With this approximation, 
equation (6.10) leads to the master equation 



Api 
'2An, 



■AAi 



(6.12) 



In other words, the frequency value that maximizes the 
.F-statistic for a given AA, does not correspond to the 
intrinsic source frequency /(g) but is instead given by a 
linear combination of the AAj. 

Let us rewrite equation (6.10) more explicitly. As 
shown in [7], can be written in terms of the ampli- 
tude modulation functions a{t) and b{t) as 



B\Fa\'' + A\Fi\' -2Cn{FaF*) 



TcohSn{f{0)) 



D 



where A, B, C, and D are constants and 



F„, = 



Tcoh/2 
Tcoh/2 



(6.13) 

x{t)a{t)e-'^^''^'^<'Ut, (6.14) 



/-'coh/'' 
x{t)b{t)e-'^^*'f'^''Ut. (6.15) 
-Tcoh/2 



Since we are interested in calculating the frequency drift 
and not the amplitude, the variation in the phase is more 
important than the amplitude modulation. Thus, the 
factors of a{t) and b{t) can be taken to be constant in the 
above equation; see [30, 31] for a discussion of the validity 
of this approximation. Thus, maximizing is equivalent 
to maximizing |X(/)|^ where X{f) is the demodulated 
Fourier transform (DeFT) defined as 

X{f) = J x{t)e-'^^''^'^^Ut. (6.16) 

With this approximation, the master equation is ob- 
tained by solving 



5|X(/,Ad;/(o),A)|2 



df 



= 0. 



(6.17) 



The details of the calculation are given in appendix A. 
The result is: 



f{t)-Fo{t) = at)-{n-nd) 



(6.18) 



where 



J^o(t) = /(o) + E^(^*)'' (6-19) 



k=l 
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and 



v(i) 



k\ 



(6.20) 



/ fd(k) t^,-,k-i \ r{t)-r{to) 

The A/(fe)'s are the residual spin-down parameters: 
^f(k) = f{k) ~ fd(k) I s-iid r(i) is the position of the de- 
tector in the SSB frame. As expected, if AA = so that 
n = nd and f(k) = fd{k), then f{t) = /(o). Further- 
more, it is clear that this master equation is qualitatively 
similar to equation (2.4) except for a constant frequency 
offset C • nrf. Thus, many of the methods obtained for the 
non-demodulated case will still be valid. 



C. Implementation details 

As mentioned above, the master equations (2.4) and 
(6.18) are qualitatively similar except for a constant fre- 
quency offset. Thus, many of the earlier results are still 
valid with some minor modifications which we now ex- 
plain. 

Resolution in parameter space: The formula for the 
resolution in space is the same as given in equation 
(4.16). However, since we can make T^oh much larger 
than before, the resolution can be made much more finer. 
Thus, for the first spin-down parameter, instead of equa- 
tion (4.17), we would have 



^/(i) = (3.7 X 10-"Hz/s) X 



365days Iday 



(6.21) 



Furthermore, the restriction due to the length of T^^i, (see 
equation (4.21)) is no longer an issue. 

As for the sky-positions, using the approximation given 
in equation (6.24), the estimate of the resolution in the 
sky proceeds in the same way as the derivation of equa- 
tion (4.12). The results of equations (4.15) and (4.16) 
are still valid, the only change being that T^^i, is now of 
the order of a day. Therefore, we rewrite equations (4.15) 
and (4.16) as: 



= -4 = 



c5J_ 
^ f 



1.0 X 10"^rad x 



/lday\ / 500 Hz 



(6.22) 

10-4^ 



v/c 



and 




(6.23) 



, r ,^2 , /^lday\ = /500Hz\ = 
= 1.5xl0-radx(^j ( 



f J \ 



10" 



The sky resolution obtained from equation (6.23) 
is therefore about 5 times better than in the non- 
demodulated case obtained from equation (4.15). 

Sky-patch size: Unlike in the non-demodulated case, 
since we are removing the frequency modulation of the 
signal beforehand, there is now, except for computational 
constraints, no restriction at all on the coherent integra- 
tion time T^oh- Typically, T^^^^ will be taken to be of 
the order of a day. However, the price we pay for this 
is that the demodulation is not valid for arbitrarily large 
patches. The patch size is determined by the largest frac- 
tional loss of sensitivity (e.g., the value) we are willing 
to tolerate from a true signal with certain mismatch pa- 
rameters A^. 

If we have demodulated for a direction in the sky, 
how different can n be from so that the loss in the 
signal power does not become unacceptably large? In 
order to answer this question, we would have to analyze 
the parameter space metric defined in terms of the mis- 
match [28] . The analysis of the metric will be presented 
elsewhere, but in this paper we shall just use a conserva- 
tive estimate for the size of the sky-patch. 

To estimate the size of the sky-patch, first note that 
the quantity ( appearing in equation (6.18) is, to a very 
good approximation, given by 



at) « m 



(6.24) 



where, as before, v is the velocity of the detector in the 
SSB frame. The velocity v(i) is the sum of two compo- 
nents, the velocity Vj, due to the yearly motion around 
the sun and the velocity due to the rotation of Earth 
around its axis: v = -|- v^. For reference, for the GEO 
detector, the magnitude Vy is about 10^ times larger than 
Vd. The estimate of the sky-patch size proceeds roughly 
like the estimate of the pixel size in equation (6.23) ex- 
cept for one difference. If we take the coherent integra- 
tion time to be roughly of the order of less than a 
day, say a third of a day, then the relevant velocity is v^. 
Thus, the sky-patch size h is roughly given by 



h: 



c_5f_ ^ c 



(6.25) 



Since Vd is roughly 100 times smaller than v, 
h w lOO(50)„i„. Thus, a typical sky-patch consists 
of about 100 pixels on a side. It should be emphasized 
that this is only an educated guess and is not likely to 
be valid for larger T^ah- 

Validity of the look up tables: Again using the ap- 
proximation given in equation (6.24), the number of fre- 
quency bins for which the LUT is valid can be estimated 
in a similar way as in the non-demodulated case. The 
master equation is 



A/ := / - Fo = /- 



nd) 



(6.26) 
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Rewrite the equation as 
1 1 
7 



V 

— -COS( 

A/ c 



1 V ■ 

A/ c 



(6.27) 



Keeping A/ fixed and differentiating with respect to / 
leads to 



dl 



1 V 
Af~c 



sin (p d4> 



(6.28) 



As before, define k and r by d/ = k5/ and dcj) = r (<5<^)^,„. 

Substituting these definitions in the above equation 
yields 



A/ 



■ sm ( 



(6.29) 



There are now two cases to look at, namely when is 
close to 7r/2 or when it is close to (or tt). First the 
easy case when (j) ^ 7r/2. Here the width of the annuli is 
roughly the same as the pixel size: 5(p ^ {S(j))^.^. Thus, 
if h is the length of a side of the sky-patch (assumed to 
be square) then the number of annuli in the sky patch is 
h/5(j) which means A/ = 5f ■ {h/5<j)). Substituting this 
in equation (6.29) and also setting <p = ■jt/2 finally leads 
to the result 



K = kq :- 



hv/c 



(6.30) 



Now turn to the large annulus case. The annulus size is 
given by Sep ~ {{5^)^,J and again Af = Sf ■ {h/Scj)). As 
for the numerator of equation (6.30), take the smallest 
value of sint/), i.e. when is no bigger than a pixel so 
that sin^ ^ {S(j))^.^. Substituting these estimates leads 
to 



(- 



Jo V 



(6.31) 



From equations (6.24) and (6.30) we see that typically, 
K for the thick annulus case is about 100 times smaller 
than for the thin annulus case. 



D. Statistics 

This subsection describes the statistics of the Hough 
map, the JF-statistic, the optimal thresholds and the sen- 
sitivity. The discussion closely parallels that of section 
V; here we simply point out some of the differences that 
arise when the ^-statistic is considered instead of the 
normalized power. 

Just as the distribution of the normalized pk power in 
subsection V A turned out to be related to the distri- 
bution with 2 degrees of freedom, one might intuitively 
expect that the distribution of should also be related 
to a distribution. However, since J-^ is constructed 
from the four filters given in equation (6.7), it turns out 



that the distribution of 2jF is a non-central distribu- 
tion with four degrees of freedom. As before, we shall 
denote the non-centrality parameter by A, and it turns 
out to be 



A = {h\h) 



(6.32) 



where the inner product (-j-) has been defined in equa- 
tion (6.2). A word of caution: while we use the same 
symbol for the non-centrality parameter as in the non- 
demodulated case, this definition is different from that of 
equation (5.3). 

Thus, the distribution of !F is 



p(^|A) = 2x'(2.F|4,A) 

1/2 



(6.33) 



= (?) /.(V^)exp(-^-i) 

where Ii is the modified Bessel function of the first order. 
In the absence of a signal, this reduces to 



p{J^\ 0) = JTe" 



(6.34) 



We select frequency bins by setting a threshold {F^i, on 
the value of the .^-'-statistic in that frequency bin. Given 
J\h, the probabilities for false alarm, false detection and 
detection are defined analogous to equations (5.6), (5.7) 
and (5.8): 



"(J^th) 

/3(:^th|A) 

^(^th|A) 



p(jc-| 0) dJ" 
(l+^,,)e"^*N 

Jo 



(6.35) 
(6.36) 

(6.37) 



The relation between a and Tt^, is different from the re- 
lation a = e"^"" in the non-demodulated case. For small 
signals, T]{Ttb\\) can be expanded as 



\ •r-2 

■n = a+ ^e-^"- + 0(A2) . 



(6.38) 



Once again we will approximate this distribution by a 

binomial. In fact, we expect the binomial approximation 
to be better in this as compared to the non-demodulated 
search because, typically, T^oh will now be larger and thus 
the signal will see a more representative 'average' of the 
detector antenna pattern. Finally, the expressions for the 
false alarm and false dismissal probabilities in the Hough 
plane are the same as in equations (5.13) and (5.14) but 
again with the A's and 77's as above. 

With the above definitions at hand, we are now ready 
to optimize the thresholds J^,h and rith using the proce- 
dure described in subsection VB. The differences from 
that subsection are simply in the dependence of a on T^-i^ 
and of r] on a. The solution for rith obtained by inverting 
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FIG. 15: Graph of I3h as a function of q = (1 +.7-th)e~'^"' for 
A = 0.2 (upper curve) and A = 0.4 (lower curve). Both curves 
correspond to an = 0.01 and A'^ = 1000. The minimum of Ph 
occurs roughly at a = 0.26 which corresponds to J^th = 2.6. 
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FIG. 16; Graph of the minimum of /3h as a function of A for 
N = 1000 and 2000. Both curves correspond to an = 0.01. 



the equation anirith, ce, N) given in figure 10 and the an- 
alytic approximation of equation (5.21) are unchanged. 
The graph of (5h as a function of a is however, now differ- 
ent. The result is shown in figure 15. The optimal value 
for the threshold turns out to be = 2.6 corresponding 
to a false alarm rate of a* = 0.26. The minimum value 
of i3h achieved by these thresholds is plotted in figure 16 
as a function of A. 

Finally, let us calculate the sensitivity of the search 
and obtain the analog of equation (5.35). The starting 
point is again equation (5.26) but now rj is related to a 
by equation (6.38) and a is related to J^th by equation 
(6.35). Then, ignoring terms of O(A^) we get the linear 
approximation for Ph- 



Ph = ^erfc 



-erfc-i(2a5,) + iee-^-(^:j2A 



where, as before, Q is given by equation (5.28). Solving 
for A in the large N limit leads to 



45 



2a*(l - a*) _ 12.73 



N 



(6.40) 



H 



where, to obtain numerical values, we have taken a 
0.01, = 0.1. Now average over the parameters 
{t, ip, a, 6) and obtain 



ho = 



8.92 



(6.41) 



In the second step we have assumed T^^, = NT^^^, which 
is valid only if there are no gaps in the data. 

As expected, equation (6.41) is identical to equation 
(5.35) except for a slightly different numerical factor. 
Thus for comparable values of Tooh and A^, the two ver- 
sions of the Hough transform search are very similar in 
sensitivity but the search with demodulated data does 
not have any restriction on T^„i, and will thus lead to a 
much greater sensitivity, though over a smaller region in 
parameter space. Thus, if we estimate the astrophysical 
range of the search as in equation (5.37), we obtain: 



TV 

d= 15.4kpcx( — 



/ fr 



VlO-e; Vldayy 



1038kg-mV V500Hz 
lO-^^Hz-i 



(6.42) 



(6.39) 



Here we have taken a coherent integration time of Iday 
and a total observation time of lyr as the reference val- 
ues. 



VII. CONCLUSIONS 

Let us summarize the main ideas and results presented 
in this work. Since it is not feasible to perform large 
parameter space searches using the matched filter with 
presently available computing power, we began by em- 
phasizing the need for hierarchical searches demonstrat- 
ing the need for an incoherent and computationally in- 
expensive search method. The Hough transform is an 
example of such a method. It looks for patterns in the 
frequency time plane by constructing a histogram in pa- 
rameter space based on the consistency of observations 
in the time-frequency plane with an underlying model 
describing the pattern. We have given a general descrip- 
tion of the Hough transform and shown its relevance for 
pulsar searches. 

We have presented two versions of the Hough trans- 
form search. The first version takes simple Fourier trans- 
forms as input data. This restricts the time baseline of 
the different segments but it allows us to search over a 
large sky-patch. The second version takes input data 
which has been demodulated to remove the effects of 
Earth's motion and the spin-down of the star; this is 
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achieved by using the jF-statistic. Wc have presented 
some technical details for both flavors of the search. In 
particular, we show how to solve the master equation in 
the two cases and how the use of look-up tables can lead 
to a large saving in computational cost. 

We have also analyzed the statistics for both cases and 
we saw that we need to choose two thresholds: the thresh- 
old Pth or J-^^^ on the coherent statistic used in the two 
cases, and the threshold on the number count in the 
Hough maps. These thresholds have been chosen in such 
a way that we get the lowest possible false dismissal rate 
for a given choice of the false alarm rate. We also estimate 
the sensitivity of the two flavors of the Hough transform 
and we find that for the same value of T^oh and N, both 
variations have comparable sensitivity, which improves 
as N~^/'^T^J^^'^ , as would be expected for an incoherent 
method that builds on coherent sub-steps. When com- 
pared to the sensitivity that a fully coherent search in a 
very large parameter space would have for the same total 
observation time T„t,j., the Hough methods arc worse by 
roughly a factor of TV^/''. Considering that the Hough 
transform can be expected to run very much faster than 
any coherent method, it should therefore be able to sur- 
vey much larger volumes of space than coherent methods, 
despite its poorer sensitivity in any single direction. This 
is therefore a potentially very important method for con- 
ducting large-scale gravitational wave pulsar surveys. 
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APPENDIX A: THE MASTER EQUATION FOR 
THE SEARCH WITH DEMODULATED DATA 

Here we detail the steps leading from equation (6.17) 
to the master equation (6.18). 

Ignoring the amplitude modulation and a possible con- 
stant phase, the signal from a pulsar with parameters 
(/(O); A) = (/(o); n, {/(„)}) would be: 



and 



/i(t,/(o),A) = e**(*'^(»)'^"), 



(Al) 



where 
$(i;/(0),A)=27r 



/(o)At„ + 



(A2) 



At„ = i3sb(i, n) - i3,b(io, n) . 



(A3) 



Here to is the time in the detector frame to which the 
frequency and spin-down parameters refer to and t^^i, is 
time in the SSB frame. Neglecting higher order relativis- 
tic effects, the detector time t is related to t^^b by 



issb(i,n) = t + 



r{t) ■ n 



(A4) 



where r{t) is the detector position in the SSB frame. 

The DeFT of the pulsar signal (Al) with respect to 
the demodulation parameters (/, A^) is: 



X{f) = 



gi[a.(t,/(o),X)-*(t,/,X<i)]^^_ 



(A5) 



2 con 



Without any loss of generality, we have taken the coher- 
ent time interval to be centered around t = so that the 
integral is from -T,„i,/2 to T,^i,/2. 

Our goal is to determine an analytical expression for 

the value of / that maximizes the power P{f) = |X(/)p 
in terms of /(o), A^ and AA. To do this, first expand 
A$(i) := $(i, /(o). A) - $(t, /, Ad) in powers of A/(fc) := 
f{k) - fd{k) and An := n - n^, keeping only terms up to 
linear order: 



^*<"-{/,.,-/)A.„.+i:7^(A'~r' 



27r 



/(0)+E4^(Ainj'^ 



fc=l 



Ar , 
An 



(A6) 



where Ar = v{t) — r(to). 

Now Taylor expand A$(f) about a fiducial time ti in 
the interval —T^„^/2 <t\< r„oh/2 again retaining terms 
only up to linear order [36] 



*(/) - /: 



e'^^'^'^dt . (A7) 



2 coh 



P(/) does not depend on A$(ti), and its maximum is 
reached for the value of / that satisfies 



9A$ 



(A8) 



f=ti 
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Differentiating (A6) with respect to t, we get r(io)- Setting the right hand side of this equation to zero 

and dropping higher order terms leads to 



1 aA$ 



27r dt 



(A9) 



\ fe=i / where Fq is as given in equation (6.19). All the dcpcn- 

fe i\ / '^(^i) \ ^i"! deuce on the residual spin-down parameters appears only 

(^ti) )("'■"' ' '^'^ ) ' ^^'^ definition of Fq and, after replacing the arbitrary 

fe=i ^ / ^ ^ time ti by we get equation (6.18). 

where Ati = issb(ii,nd) - Cb(io,nd) and Ari = r{ti) - 
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